function qs = initialize_homogeneous_subsidies(mu,B,A,P,rho,lambda_rnd,lambda_tech)

%% Computation of the M matrix.
n = length(A);
M = zeros(n,n);
q_bar = zeros(n,1);
r = zeros(n,1);
u = ones(n,1);
M = inv(eye(n)+rho*B-lambda_rnd*A-lambda_tech*P);
q_bar = M*mu; % Vector of weighted Bonacich centralities.
r = M*(eye(n)+lambda_rnd*A+lambda_tech*P)*u;

%% Compute the homogeneous subsidy.
numerator = 0;
denominator = 0;
for i=1:n
    dummy1 = 0;
    dummy2 = 0;
    for j=1:n
        dummy1 = dummy1 + B(i,j)*(q_bar(i)*r(j)+q_bar(j)*r(i));
        dummy2= dummy2 + B(i,j)*r(j);
    end
    numerator = numerator + q_bar(i)*(2*r(i)-1) + rho/2*dummy1;
    denominator = denominator + 1 + r(i)*(2 - 2*r(i)-rho*dummy2);
end
subsidy = numerator/denominator;

%% Compute output levels with the subsidy.
q_sub = q_bar + subsidy*r;

%% Check if output levels are all non-negative. If not, determine corner solution.
ind = find(q_sub<0);
if(~isempty(ind))
    q_sub = output_corner(mu,B,A,P,rho,lambda_rnd,lambda_tech,subsidy*ones(n,1));
end

%% Check if zero subsidies generate higher welfare.
q_bar = output_corner(mu,B,A,P,rho,lambda_rnd,lambda_tech,zeros(n,1));
W = q_bar'*q_bar + rho/2*q_bar'*B*q_bar;
W_sub = q_sub'*q_sub + rho/2*q_sub'*B*q_sub - subsidy*sum(q_sub) - n/2*subsidy^2;
if(W>W_sub)
    q_sub = q_bar;
    subsidy = 0;
end

qs = [q_sub; subsidy];
